#setwd('/afs/inf.ed.ac.uk/user/s17/s1725186/Documents/PhD-Models/FirstPUModel/RMarkdowns')

library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(plotlyutils)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(GGally)
library(Rtsne)
library(ClusterR)
library(DESeq2)
library(expss)
library(knitr)

Load preprocessed dataset (preprocessing code in 19_10_14_data_preprocessing.Rmd)

# Gandal dataset
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
DE_info = DE_info %>% data.frame

# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID'=as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal'=1)

# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_01-03-2020_w_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]

# Update DE_info with SFARI and Neuronal information
DE_info = DE_info %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>% 
  mutate(`gene-score`=ifelse(is.na(`gene-score`), 'Others', `gene-score`)) %>%
  distinct(ID, .keep_all = TRUE) %>% left_join(GO_neuronal, by='ID') %>%
  mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
  mutate(gene.score=ifelse(`gene-score`=='Others' & Neuronal==1, 'Neuronal', `gene-score`), significant=padj<0.05 & !is.na(padj)) %>%
  mutate(Group = factor(ifelse(gene.score %in% c('Neuronal','Others'), gene.score, 'SFARI'), levels = c('SFARI', 'Neuronal', 'Others')))


SFARI_colour_hue = function(r) {
  pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}

rm(GO_annotations)

SFARI Gene list

cat(paste0('There are ', length(unique(SFARI_genes$`gene-symbol`)), ' genes with a SFARI score'))
## There are 912 genes with a SFARI score

There are 912 genes with a SFARI score. but to map them to gene expression mapa we had to map the gene names to their corresponding ensembl IDs

Mapping SFARI Gene names to Ensembl IDs

cat(paste0('There are ', nrow(SFARI_genes), ' Ensembl IDs corresponding to the ',
             length(unique(SFARI_genes$`gene-symbol`)),' genes in the SFARI Gene dataset'))
## There are 1116 Ensembl IDs corresponding to the 912 genes in the SFARI Gene dataset



  • Since a gene can have more than one ensembl ID, there were some one-to-many mappings between a gene name and ensembl IDs, so that’s why we ended up with 1090 rows in the SFARI_genes dataset.

  • The details about how the genes were annotated with their Ensembl IDs can be found in 20_02_06_get_ensembl_ids.html

cat(paste0('There are ', sum(is.na(SFARI_genes$`gene-score`)) ,
             ' genes in the SFARI list without a score, of which ',
             sum(is.na(SFARI_genes$`gene-score`) & SFARI_genes$syndromic==0),
             ' don\'t have syndromic tag either'))
## There are 87 genes in the SFARI list without a score, of which 0 don't have syndromic tag either

Exploratory Analysis

cat(paste0('There are ', sum(SFARI_genes$ID %in% rownames(datExpr)), ' SFARI Genes in the expression dataset (~',
             round(100*mean(SFARI_genes$ID %in% rownames(datExpr))),'%)'))
## There are 864 SFARI Genes in the expression dataset (~77%)
cat(paste0('Of these, only ', sum(DE_info$`gene-score`!='Others'), ' have an assigned score'))
## Of these, only 789 have an assigned score

From now on, we’re only going to focus on the 789 genes with a score

Gene count by SFARI score:

table_info = DE_info %>% apply_labels(`gene-score` = 'SFARI Gene Score', syndromic = 'Syndromic Tag',
                                      Neuronal = 'Neuronal Function', gene.score = 'Gene Score') %>%
             mutate(syndromic = as.logical(syndromic), Neuronal = as.logical(Neuronal))

cro(table_info$`gene-score`)
 #Total 
 SFARI Gene Score 
   1  144
   2  205
   3  440
   Others  15358
   #Total cases  16147

Gene count by Syndromic tag:

cro(table_info$syndromic)
 #Total 
 Syndromic Tag 
   FALSE  748
   TRUE  116
   #Total cases  864


GO Neuronal annotations:

cat(glue(sum(GO_neuronal$ID %in% rownames(datExpr)), ' genes have neuronal-related annotations'))
## 1094 genes have neuronal-related annotations
cat(glue(sum(SFARI_genes$ID %in% GO_neuronal$ID),' of these genes have a SFARI score'))
## 179 of these genes have a SFARI score
cro(table_info$gene.score[DE_info$`gene-score` %in% c('1','2','3','4','5','6')],
    list(table_info$Neuronal[DE_info$`gene-score` %in% c('1','2','3','4','5','6')], total()))
 Neuronal Function     #Total 
 FALSE   TRUE   
 Gene Score 
   1  105 39   144
   2  161 44   205
   3  365 75   440
   #Total cases  631 158   789



All SFARI scores together


Gene Expression


Larger mean expression than the other two groups, smaller SD than Neuronal genes

plot_data = data.frame('ID'=rownames(datExpr), 'MeanExpr'=rowMeans(datExpr), 'SDExpr'=apply(datExpr,1,sd)) %>% 
            left_join(DE_info, by='ID') %>% mutate(Group = factor(ifelse(gene.score %in% c('Neuronal','Others'), gene.score, 'SFARI'), 
                                                                  levels = c('SFARI', 'Neuronal', 'Others')))

p1 = ggplotly(plot_data %>% ggplot(aes(Group, MeanExpr, fill=Group)) + geom_boxplot() + 
              scale_fill_manual(values=c('#00A4F7', SFARI_colour_hue(r=c(8,7)))) + theme_minimal() +
              theme(legend.position='none'))

p2 = ggplotly(plot_data %>% ggplot(aes(Group, SDExpr, fill=Group)) + geom_boxplot() + 
              scale_fill_manual(values=c('#00A4F7', SFARI_colour_hue(r=c(8,7)))) + theme_minimal() +
              ggtitle('Mean Expression (left) and SD (right) comparison between SFARI Genes and the rest of the dataset') + 
              theme(legend.position='none'))

subplot(p1, p2, nrows=1)
rm(plot_data, p1, p2)

Log Fold Change

Proportion of over- and under-expressed genes is very similar between groups: approximately half

DE_info %>% mutate(direction = ifelse(log2FoldChange>0, 'over-expressed', 'under-expressed')) %>% group_by(Group, direction) %>% tally(name = 'over_expr') %>% 
            filter(direction == 'over-expressed') %>% ungroup %>% left_join(DE_info %>% group_by(Group) %>% tally(name = 'Total'), by = 'Group') %>% ungroup %>%
            mutate('prop_over_expr' = round(over_expr/Total,3)) %>% dplyr::select(-direction) %>% kable
Group over_expr Total prop_over_expr
SFARI 380 789 0.482
Neuronal 461 936 0.493
Others 7419 14422 0.514

Smaller lFC magnitude than Neuronal genes and similar but slightly lower than the rest of the genes

ggplotly(DE_info %>% ggplot(aes(x=Group, y=abs(log2FoldChange), fill=Group)) + geom_boxplot() + 
         scale_fill_manual(values=c('#00A4F7', SFARI_colour_hue(r=c(8,7)))) + ylab('logFoldChange Magnitude') + xlab('Group') + theme_minimal() +
         theme(legend.position='none'))
  • SFARI Genes, as a group, have less genes with high (positive) lFC than the rest of the genes in the dataset

  • Perhaps the opposite is true for the genes with the highest (negative) lFC, although this pattern is weaker

plot_data = DE_info %>% mutate(Group = ifelse(gene.score %in% c('Neuronal', 'Others'), gene.score, 'SFARI')) %>% dplyr::select(Group, log2FoldChange) %>%
            mutate(Group = factor(Group, levels = c('Others', 'Neuronal', 'SFARI')),
                   quant = cut(log2FoldChange, breaks = quantile(log2FoldChange, probs = seq(0,1,0.05)) %>% as.vector, labels = FALSE)) %>% 
            filter(Group == 'SFARI') %>% group_by(quant) %>% tally %>% ungroup

ggplotly(plot_data %>% ggplot(aes(quant, n)) + geom_bar(stat = 'identity', fill = '#00A4F7') + geom_smooth(color = 'gray', alpha = 0.1) + 
         xlab('Log Fold Change Quantiles') + ylab('Number of SFARI Genes') + ggtitle('Number of SFARI Genes for lFC Quantiles') + theme_minimal())
cat('LFC values for each quantile:')
## LFC values for each quantile:
quants = data.frame('Quantile' = 1:20, 'LFC Range' = cut(DE_info$log2FoldChange, breaks = quantile(DE_info$log2FoldChange, probs = seq(0,1,0.05)) %>% as.vector) %>%
         table %>% names)

quants %>% kable
Quantile LFC.Range
1 (-2.04,-0.229]
2 (-0.229,-0.164]
3 (-0.164,-0.125]
4 (-0.125,-0.0988]
5 (-0.0988,-0.0777]
6 (-0.0777,-0.0598]
7 (-0.0598,-0.0427]
8 (-0.0427,-0.0261]
9 (-0.0261,-0.0114]
10 (-0.0114,0.00376]
11 (0.00376,0.0188]
12 (0.0188,0.0352]
13 (0.0352,0.052]
14 (0.052,0.0731]
15 (0.0731,0.097]
16 (0.097,0.129]
17 (0.129,0.171]
18 (0.171,0.238]
19 (0.238,0.373]
20 (0.373,1.96]
rm(quants)


DEA

Lower proportion of DE genes than genes with Neuronal annotation but higher than the rest of the genes

DE_info %>% group_by(Group, significant) %>% tally(name = 'DE') %>% filter(significant) %>% ungroup %>%
            left_join(DE_info %>% group_by(Group) %>% tally(name = 'Total'), by = 'Group') %>% ungroup %>% mutate('prop_DE' = round(DE/Total,2)) %>% 
            dplyr::select(-significant) %>% kable
Group DE Total prop_DE
SFARI 244 789 0.31
Neuronal 347 936 0.37
Others 3908 14422 0.27

SFARI Genes have consistently a lower percentage of DE genes than the Genes with Neuronal annotations and a very similar level to the rest of the genes

The SFARI Genes have less DE genes than when using the original SFARI genes (makes sense since we are losing all genes with SFARI scores 5 and 6, which were the ones with the largest lFC)

lfc_list = seq(1, 1.2, 0.01)

all_counts = data.frame('group'='All', 'n'=as.character(nrow(DE_info)))
Others_counts = data.frame('group'='Others', n=as.character(sum(DE_info$Group == 'Others')))
Neuronal_counts = data.frame('group'='Neuronal', n=as.character(sum(DE_info$Neuronal)))
lfc_counts_all = DE_info %>% filter(Group == 'SFARI') %>% tally %>%
                 mutate('group'='SFARI', 'n'=as.character(n)) %>%
                 dplyr::select(group, n) %>%
                 bind_rows(Neuronal_counts, Others_counts, all_counts) %>%
                 mutate('lfc'=-1) %>%  dplyr::select(lfc, group, n)

for(lfc in lfc_list){
  
  # Recalculate DE_info with the new threshold (p-values change)
  DE_genes = results(dds, lfcThreshold=log2(lfc), altHypothesis='greaterAbs') %>% data.frame %>% mutate('ID'=rownames(.)) %>% 
             left_join(DE_info %>% dplyr::select(ID, Neuronal, gene.score, Group), by = 'ID') %>% filter(padj<0.05 & abs(log2FoldChange)>log2(lfc))

  
  # Calculate counts by groups
  all_counts = data.frame('group'='All', 'n'=as.character(nrow(DE_genes)))
  Others_counts = data.frame('group'='Others', n=as.character(sum(DE_genes$Group == 'Others')))
  Neuronal_counts = data.frame('group'='Neuronal', n=as.character(sum(DE_genes$Neuronal)))
  lfc_counts = DE_genes %>% filter(Group == 'SFARI') %>% tally %>%
               mutate('group'='SFARI', 'n'=as.character(n)) %>%
               bind_rows(Neuronal_counts, Others_counts, all_counts) %>%
               mutate('lfc'=lfc) %>% dplyr::select(lfc, group, n)
  
  
  # Update lfc_counts_all
  lfc_counts_all = lfc_counts_all %>% bind_rows(lfc_counts)
}

# Add missing entries with 0s
lfc_counts_all = expand.grid('group'=unique(lfc_counts_all$group), 'lfc'=unique(lfc_counts_all$lfc)) %>% 
  left_join(lfc_counts_all, by=c('group','lfc')) %>% replace(is.na(.), 0)

# Calculate percentage of each group remaining
tot_counts = DE_info %>% filter(Group == 'SFARI') %>% tally() %>%
             mutate('group'='SFARI', 'tot'=n) %>% dplyr::select(group, tot) %>%
             bind_rows(data.frame('group'='Neuronal', 'tot'=sum(DE_info$Neuronal)),
                       data.frame('group' = 'Others', 'tot' = sum(DE_info$Group == 'Others')),
                       data.frame('group'='All', 'tot'=nrow(DE_info)))

lfc_counts_all = lfc_counts_all %>% filter(lfc!=-1) %>% #, group!='Others') %>% 
                 left_join(tot_counts, by='group') %>% mutate('perc'=round(100*as.numeric(n)/tot,2))


# Plot change of number of genes
ggplotly(lfc_counts_all %>% filter(group != 'All') %>% mutate(group = factor(group, levels = c('SFARI', 'Neuronal', 'Others'))) %>%
         ggplot(aes(lfc, perc, color=group)) + geom_point(aes(id=n)) + geom_line() +
         scale_color_manual(values=c('#00A4F7', SFARI_colour_hue(r=c(8,7)))) + ylab('% of Differentially Expressed Genes') +  xlab('Fold Change') +
         ggtitle('Effect of filtering thresholds in SFARI Genes') + theme_minimal())



Grouping Genes by SFARI Gene Score


Gene Expression




Normalised data

  • The higher the SFARI score, the higher the mean expression of the gene: This pattern is quite strong and it doesn’t have any biological interpretation, so it’s probably bias in the SFARI score assignment

  • The higher the SFARI score, the lower the standard deviation: This pattern is not as strong, but it is weird because the data was originally heteroscedastic with a positive relation between mean and variance, so the fact that the relation now seems to have reversed could mean that the vst normalisation ended up affecting the highly expressed genes more than it should have when trying to correct their higher variance

plot_data = data.frame('ID'=rownames(datExpr), 'MeanExpr'=rowMeans(datExpr), 'SDExpr'=apply(datExpr,1,sd)) %>% 
            left_join(DE_info, by='ID')

p1 = ggplotly(plot_data %>% ggplot(aes(gene.score, MeanExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + theme_minimal() +
              theme(legend.position='none'))

p2 = ggplotly(plot_data %>% ggplot(aes(gene.score, SDExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + theme_minimal() +
              ggtitle('Mean Expression (left) and SD (right) by SFARI score') + 
              theme(legend.position='none'))

subplot(p1, p2, nrows=1)
rm(plot_data, p1, p2)


Raw data

Just to corroborate that the relation between sd and SFARI score used to be in the opposite direction before the normalisation: The higher the SFARI score the higher the mean expression and the higher the standard deviation

*There are a lot of outliers, but the plot is interactive so you can zoom in

# Save preprocessed results
datExpr_prep = datExpr
datMeta_prep = datMeta
DE_info_prep = DE_info

load('./../Data/filtered_raw_data.RData')

plot_data = data.frame('ID'=rownames(datExpr), 'MeanExpr'=rowMeans(datExpr), 'SDExpr'=apply(datExpr,1,sd)) %>% 
            left_join(DE_info, by='ID')

p1 = ggplotly(plot_data %>% ggplot(aes(gene.score, MeanExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + theme_minimal() +
              theme(legend.position='none'))

p2 = ggplotly(plot_data %>% ggplot(aes(gene.score, SDExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + theme_minimal() +
              ggtitle('Mean Expression (left) and SD (right) by SFARI score') + 
              theme(legend.position='none'))

subplot(p1, p2, nrows=1)
rm(plot_data, p1, p2)

Return to normalised version of the data

# Save preprocessed results
datExpr = datExpr_prep
datMeta = datMeta_prep
DE_info = DE_info_prep

rm(datExpr_prep, datMeta_prep, DE_info_prep)


Log Fold Change

There seems to be a negative relation between SFARI score and log fold change when it would be expected to be either positively correlated or independent from each other (this last one because there are other factors that determine if a gene is releated to Autism apart from differences in gene expression)

Wikipedia mentions the likely explanation for this: “A disadvantage and serious risk of using fold change in this setting is that it is biased and may misclassify differentially expressed genes with large differences (B − A) but small ratios (B/A), leading to poor identification of changes at high expression levels”.

Based on this, since we saw there is a strong relation between SFARI score and mean expression, the bias in log fold change affects mainly genes with high SFARI scores, which would be the ones we are most interested in.

On top of this, I believe this effect is made more extreme by the pattern found in the previous plots, since the higher expressed genes were the most affected by the normalisation transformation, ending up with a smaller variance than the rest of the data, which is related to smaller ratios. (This is a constant problem independently of the normalisation function used).

ggplotly(DE_info %>% ggplot(aes(x=gene.score, y=abs(log2FoldChange), fill=gene.score)) + 
         geom_boxplot() + scale_fill_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + xlab('SFARI Gene scores') + ylab('log Fold Change Magnitude') +
         theme_minimal() + theme(legend.position='none'))


Effects of modifying filtering threshold by SFARI score

  • The % of DE Genes is very similar for all 3 scores and for genes with no annotation

  • All SFARI scores except are more affected by the filtering than the genes with Neuronal annotations

  • Using the null hypothesis \(H_0: lfc=0\), 244/789 SFARI genes are statistically significant (31%)

lfc_list = seq(1, 1.2, 0.01)

all_counts = data.frame('group'='All', 'n'=as.character(nrow(DE_info)))
Neuronal_counts = data.frame('group'='Neuronal', n=as.character(sum(DE_info$Neuronal)))
lfc_counts_all = DE_info %>% group_by(`gene-score`) %>% tally %>%
                 mutate('group'=as.factor(`gene-score`), 'n'=as.character(n)) %>%
                 dplyr::select(group, n) %>%
                 bind_rows(Neuronal_counts, all_counts) %>%
                 mutate('lfc'=-1) %>%  dplyr::select(lfc, group, n)

for(lfc in lfc_list){
  
  # Recalculate DE_info with the new threshold (p-values change)
  DE_genes = results(dds, lfcThreshold=log2(lfc), altHypothesis='greaterAbs') %>% data.frame
  
  DE_genes = DE_genes %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>% 
             mutate(`gene-score`=ifelse(is.na(`gene-score`), 'Others', `gene-score`)) %>%
             distinct(ID, .keep_all = TRUE) %>% left_join(GO_neuronal, by='ID') %>%
             mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
             mutate(gene.score=ifelse(`gene-score`=='Others' & Neuronal==1, 'Neuronal', `gene-score`))
  
  DE_genes = DE_genes %>% filter(padj<0.05 & abs(log2FoldChange)>log2(lfc))

  
  # Calculate counts by groups
  all_counts = data.frame('group'='All', 'n'=as.character(nrow(DE_genes)))
  Others_counts = data.frame('group'='Others', n=as.character(sum(DE_genes$gene.score == 'Others')))
  Neuronal_counts = data.frame('group'='Neuronal', n=as.character(sum(DE_genes$Neuronal)))
  lfc_counts = DE_genes %>% group_by(`gene-score`) %>% tally %>%
               mutate('group'=`gene-score`, 'n'=as.character(n)) %>%
               bind_rows(Neuronal_counts, all_counts) %>%
               mutate('lfc'=lfc) %>% dplyr::select(lfc, group, n)
  
  
  # Update lfc_counts_all
  lfc_counts_all = lfc_counts_all %>% bind_rows(lfc_counts)
}

# Add missing entries with 0s
lfc_counts_all = expand.grid('group'=unique(lfc_counts_all$group), 'lfc'=unique(lfc_counts_all$lfc)) %>% 
  left_join(lfc_counts_all, by=c('group','lfc')) %>% replace(is.na(.), 0)

# Calculate percentage of each group remaining
tot_counts = DE_info %>% group_by(`gene-score`) %>% tally() %>% filter(`gene-score`!='Others') %>%
             mutate('group'=`gene-score`, 'tot'=n) %>% dplyr::select(group, tot) %>%
             bind_rows(data.frame('group'='Neuronal', 'tot'=sum(DE_info$Neuronal)),
                       data.frame('group' = 'Others', 'tot' = sum(DE_info$gene.score == 'Others')),
                       data.frame('group'='All', 'tot'=nrow(DE_info)))

lfc_counts_all = lfc_counts_all %>% filter(lfc!=-1) %>% #, group!='Others') %>% 
                 left_join(tot_counts, by='group') %>% mutate('perc'=round(100*as.numeric(n)/tot,2))


# Plot change of number of genes
ggplotly(lfc_counts_all %>% filter(group != 'All') %>% ggplot(aes(lfc, perc, color=group)) + geom_point(aes(id=n)) + geom_line() + 
         scale_color_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + ylab('% of Differentially Expressed Genes') +  xlab('Fold Change') + 
         ggtitle('Effect of filtering thresholds by SFARI score') + theme_minimal())
rm(lfc_list, all_counts, Neuronal_counts, lfc_counts_all, lfc, lfc_counts, lfc_counts_all, tot_counts, lfc_counts_all, Others_counts)

Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] knitr_1.28                  expss_0.10.2               
##  [3] DESeq2_1.24.0               SummarizedExperiment_1.14.1
##  [5] DelayedArray_0.10.0         BiocParallel_1.18.1        
##  [7] matrixStats_0.56.0          Biobase_2.44.0             
##  [9] GenomicRanges_1.36.1        GenomeInfoDb_1.20.0        
## [11] IRanges_2.18.3              S4Vectors_0.22.1           
## [13] BiocGenerics_0.30.0         ClusterR_1.2.1             
## [15] gtools_3.8.2                Rtsne_0.15                 
## [17] GGally_1.5.0                gridExtra_2.3              
## [19] viridis_0.5.1               viridisLite_0.3.0          
## [21] RColorBrewer_1.1-2          plotlyutils_0.0.0.9000     
## [23] plotly_4.9.2                glue_1.3.2                 
## [25] reshape2_1.4.3              forcats_0.5.0              
## [27] stringr_1.4.0               dplyr_0.8.5                
## [29] purrr_0.3.3                 readr_1.3.1                
## [31] tidyr_1.0.2                 tibble_3.0.0               
## [33] ggplot2_3.3.0               tidyverse_1.3.0            
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       ellipsis_0.3.0         htmlTable_1.13.3      
##  [4] XVector_0.24.0         base64enc_0.1-3        fs_1.4.0              
##  [7] rstudioapi_0.11        bit64_0.9-7            AnnotationDbi_1.46.1  
## [10] fansi_0.4.1            lubridate_1.7.4        xml2_1.2.5            
## [13] splines_3.6.3          geneplotter_1.62.0     Formula_1.2-3         
## [16] jsonlite_1.6.1         annotate_1.62.0        broom_0.5.5           
## [19] cluster_2.1.0          dbplyr_1.4.2           png_0.1-7             
## [22] compiler_3.6.3         httr_1.4.1             backports_1.1.5       
## [25] assertthat_0.2.1       Matrix_1.2-18          lazyeval_0.2.2        
## [28] cli_2.0.2              acepack_1.4.1          htmltools_0.4.0       
## [31] tools_3.6.3            gmp_0.5-13.6           gtable_0.3.0          
## [34] GenomeInfoDbData_1.2.1 Rcpp_1.0.4             cellranger_1.1.0      
## [37] vctrs_0.2.4            nlme_3.1-147           crosstalk_1.1.0.1     
## [40] xfun_0.12              rvest_0.3.5            lifecycle_0.2.0       
## [43] XML_3.99-0.3           zlibbioc_1.30.0        scales_1.1.0          
## [46] hms_0.5.3              yaml_2.2.1             memoise_1.1.0         
## [49] rpart_4.1-15           RSQLite_2.2.0          reshape_0.8.8         
## [52] latticeExtra_0.6-29    stringi_1.4.6          highr_0.8             
## [55] genefilter_1.66.0      checkmate_2.0.0        rlang_0.4.5           
## [58] pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14         
## [61] lattice_0.20-41        labeling_0.3           htmlwidgets_1.5.1     
## [64] bit_1.1-15.2           tidyselect_1.0.0       plyr_1.8.6            
## [67] magrittr_1.5           R6_2.4.1               generics_0.0.2        
## [70] Hmisc_4.4-0            DBI_1.1.0              mgcv_1.8-31           
## [73] pillar_1.4.3           haven_2.2.0            foreign_0.8-76        
## [76] withr_2.1.2            survival_3.1-12        RCurl_1.98-1.1        
## [79] nnet_7.3-14            modelr_0.1.6           crayon_1.3.4          
## [82] rmarkdown_2.1          jpeg_0.1-8.1           locfit_1.5-9.4        
## [85] grid_3.6.3             readxl_1.3.1           data.table_1.12.8     
## [88] blob_1.2.1             reprex_0.3.0           digest_0.6.25         
## [91] xtable_1.8-4           munsell_0.5.0